backward_differentiation_coefficients(1, 1)
backward_differentiation_coefficients(2, 2)
backward_differentiation_coefficients(3, 3)
backward_differentiation_coefficients(4, 4)

function backward_differentiation_coefficients(s, p)
    % 定义符号变量 a1, a2, ..., a_s (对应a0,a1,...,a_{s-1}) 和 bs
    syms a [1, s]    % a1对应a0, a2对应a1,...,as对应a_{s-1}
    syms bs
    syms as
    
    % 设置阶条件方程
    eqns = sym(zeros(1, p+1));  % 需要 p 个方程
    for m = 0:p  % 从0到p-1阶条件
        j = 0:s-1; % j = 0,1,...,s-1
        if m == 0
            % C0条件: sum(a_k) + 1 = 0
            eqns(1) = sum(a) + 1 == 0;
        else
            % 高阶条件: sum(j^m * a_j) + s^m = m*s^{m-1}*bs
            eqns(m+1) = sum(j.^m .* a) + s^m == m*s^(m-1)*bs;
        end
    end
    eqn = as == 1;
    % 求解方程组
    sol = solve(eqns, [a, bs]);
    sols = solve(eqn, as);
    % 提取并显示系数（转换为分数形式）
    fprintf('\nBackward Differentiation 系数 (s=%d, p=%d):\n', s, p);
    
    % 处理不同的解格式
    if s == 1
        % s=1特殊情况
        fprintf('a0 = %s\n', char(rats(double(sol.a1))));
        fprintf('a1 = %s\n', char(rats(double(sols))));
        fprintf('bs = %s\n', char(rats(double(sol.bs))));
    else
        % 一般情况
        for k = 1:s
            ak = sol.(sprintf('a%d', k));
            fprintf('a%d = %s\n', k-1, char(rats(double(ak))));
        end
        fprintf('a%d = %s\n',s, char(rats(double(sols))));
        fprintf('bs = %s\n', char(rats(double(sol.bs))));
    end
end